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I. INTRODUCTION 

Relativistic heavy ion collisions provide the means 
to study both the novel properties of the quark gluon 
plasma and the fascinating nature of how it is created 
and evolves. Unfortunately, experimental procedure is 
confined to measurements of the asymptotic momenta of 
the particles comprising the collision's debris. Addressing 
the fundamental questions concerning the properties of 
super-hadronic matter and the collision's evolution inher- 
ently depends on large-scale multistage transport models. 
Such models have improved significantly in recent years, 
and now typically combine viscous hydrodynamic treat- 
ments for the evolution of the semi-thermalized quark- 
gluon plasma (~l-7 fm/c), and microscopic hadronic sim- 
ulations to describe the dissolution and breakup of the 
produced hadrons (~7-20 fm/c). For the first fm/c of 
the collision, when the system is too far from equilibrium 
for even a viscous hydrodynamic treatment, quantitative 
modeling carries large uncertainties. If the profile and 
flow of the matter being fed into the hydrodynamic treat- 
ment could be determined phenomenologically, it would 
be invaluable in understanding how QCD saturation phe- 
nomena affect the stopping of the matter, and in deter- 
mining the mechanism and timescale for thermalization. 

The data sets from the Relativistic Heavy Ion Collider 
(RHIC) and from the heavy ion programs at the Large 



Hadron Collider (LHC) are immense. The heterogenous 
nature of the data, along with the strong interdependence 
of disparate observables with respect to basic model pa- 
rameters, makes any interpretation of the data challeng- 
ing. The phenomenology of heavy ion collisions has pro- 
gressed despite these difficulties, primarily by identifying 
the principal connections between model parameters and 
observables. For example, it is well understood that the 
shear viscosity of the quark-gluon plasma strongly affects 
the observed anisotropic flow coefficients. In an early 
analysis pQ, the viscosity is adjusted until one finds a sat- 
isfactory fit with the anisotropic flow coefficient v^- The 
shortcoming of such an approach is that several other 
unknown parameters, such as the spatial anisotropy of 
the initial state [26], also affect V2- In turn, each of these 
parameters also affects numerous other observables. Sim- 
ilar approaches with more advanced models [2HT0] have 
considered the variation of several parameters, and also 
the effects of such parameters on spectra. However, due 
to the numerical costs of running the models, these ap- 
proaches have been unable to consider the simultaneous 
variation of more than two or three parameters, or to con- 
sider a wider range of experimental observables. These 
limitations compromise both the rigor and completeness 
of the effort. 

Other fields of science face similar challenges, a no- 
table case is the extraction of cosmological parameters 
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from observations of fluctuations of the cosmic microwave 
background [TTJ [12] . Here the parameters are some of the 
most fundamental in nature, such as the densities of dark 
matter and of dark energy. To overcome the limitations 
of running the model a large number of times, a surro- 
gate model (a.k.a. an emulator) was developed to stand 
in for the true computer code. Rather than re-running 
the full cosmological evolution model during the explo- 
ration of the parameter space, one runs the full model 
at only ~ 100 — 200 points in parameter space, carefully 
chosen to best fill the overall space. A surrogate model 
was constructed that effectively interpolated from the fi- 
nite set of observations of the full model. The emulator 
was then substituted for the full model when exploring 
the parameter space . Similar ideas have been applied 
to the field of galaxy formation [13]. Here, we report 
first results for a large scale surrogate-model-based sta- 
tistical analysis of heavy-ion collision data, a small scale 
application of these ideas was discussed in [13] . 

For this first effort, only a small subset of possible data 
will be considered, that coming from 100 A GeV +100 A 
GeV Au+Au collisions at RHIC. Spectra for pions, kaons 
and protons will be considered along with the elliptic 
flow observable measured for pions, and femtoscopic 
source radii from two-pion correlations. The motivation 
for first considering soft observables is two- fold. First, 
they are the most sensitive to the model parameters re- 
lated to the bulk properties of matter, and secondly, the 
dependencies are highly intertwined. During the last two 
years, the data set for relativist ic heavy ion collisions has 
greatly expanded with the beam-energy scan at RHIC, 
and with the inaugural heavy- ion run at the LHC. The set 
will grow again in the next year as data is analyzed from 
Cu+Au and from U+U runs at RHIC. Ultimately, one 
may wish to incorporate other observables, such as dilep- 
ton emission, higher flow moments, or long-range correla- 
tions, once the theoretical treatments become more stan- 
dardized and robust. The methods described here should 
scale well with increasingly large data sets, and incorpo- 
rating additional observables into the analysis should be 
tractable. 

On the theory side, numerous parameters factor into 
models of heavy- ion collisions. Several of these param- 
eters are needed to describe the initial energy density 
and flow profiles the comprise the initial state of the hy- 
drodynamic evolution. Other parameters describe the 
bulk properties of super-hadronic matter, such as the 
equation of state and viscosities. Still other parameters 
could describe out-of-equilibrium behavior such as chem- 
ical abundances of various quark species. For tract abil- 
ity only a half dozen parameters will be considered for 
this study. Four of the parameters describe the initial 
state for the hydrodynamic module, and two describe 
the shear viscosity and its energy dependence above the 
transition temperature. The equation of state from lat- 
tice calculations Q~5] [16] will be assumed to be correct. 
In a future study, that too will be parameterized to learn 
to what extent the equation of state is constrained ex- 



perimentally. Hadronization will be assumed to produce 
a chemically-equilibrated hadronic gas when the energy 
density reaches 170 MeV. In the future, this assumption 
will also be relaxed and the away-from-equilibrium prop- 
erties of these hadrons will be parameterized. Addition- 
ally, one should expect a non-negligible bulk viscosity in 
the transition region [T7J[l8]. However, due to some nu- 
merical instabilities with handling bulk viscosity, which 
will be set to zero for this study. An advantage of sur- 
rogate model techniques is that they scale well with an 
increasing number of parameters, and the efficiency of the 
methods are not greatly diminished if several parameters 
have only marginal impact. We expect these methods to 
continue to work even if we triple the number of param- 
eters. 

Details of the model and data used for the analysis 
are provided in the next two sections. The theory of the 
model emulator is described in Sec. IV, with a test of 
the emulator against a mock data set in the subsequent 
section. Results from an analysis of the real data set are 
given in Sec. VII, while a summary and outlook comprise 
the final section. 



II. MODELING THE EVOLUTION AT RHIC 

For this study, four elements are involved in the mod- 
eling. 

1. The pre-thermal, or stopping stage. Rather than 
dynamically solving for the evolution during this 
stage, we apply a parameterized description of the 
stress-energy tensor describing the state of the col- 
lision at a time of 0.8 fm/c. Although sophisticated 
models of the initial state do exist, e.g. [T9H23] , the 
large uncertainties and the lack of theoretical con- 
sensus dissuades one from picking any individual 
model. 

2. The hydrodynamic stage lasts from 0.8 fm/c until 
the system falls below a hadronization temperature 
of 170 MeV. Viscous hydrodynamics is justified for 
a strongly interacting system that is not too far 
from equilibrium, and is especially convenient for a 
system undergoing a transition in degrees of free- 
dom, because the equations can be applicable even 
when there are no well-defined quasi-particles. 

3. Once the density has fallen to the point that the 
evolution can be modeled as binary collisions of 
hadrons, we switch to a microscopic hadronic sim- 
ulation, or cascade. The cascade is able to handle 
the loss of equilibrium between species, e.g., the 
protons and pions moving with different average 
velocities or having different kinetic temperatures. 
The cascade also handles disassociation seamlessly. 

4. Particles are correlated at small relative momen- 
tum due to interactions in the asymptotic state. 
Assuming that interactions with third bodies are 
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randomizing, one can calculate two-particle corre- 
lations given the source function, which describes 
the relative distance between two particles of simi- 
lar velocities. Taking the source function from the 
information about last collisions in the cascade, and 
convoluting with the square of the known outgo- 
ing relative two-particle wave functions, we calcu- 
late correlations, and from the correlations calcu- 
late effective Gaussian source radii which can be 
compared to those extracted from experimentally 
measured correlations functions. 



A. Parameterizing the initial state 



describes the flow profile. The first is a weight, f wn be- 
tween a wounded- nucleon profile and a saturation-based 
profile, 

e(x, y) = /wnCwnO, V) + (1 - /wnKatfc, V)- (1) 

The wounded-nucleon profile [24] and the saturation pro- 
files are based on Glauber thickness functions which 
describe the projected areal densities of the incoming 
nucelei in a plane perpendicular to the beam axis, 



T AjB (x,y) = J dz p AiB (x,y,z), 



(2) 



Rather than applying one of the competing models for 
the initial state, a parameterized form is used for the 
initial energy-density and flow profiles. Three parame- 
ters describe the initial energy-density profile and one 



where Pa,b are the baryon densities of the two nuclei 
given the impact parameters. The thickness functions 
have units of baryons per fm 2 , and the energy densities 
have the form, 
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(4) 



Here, the energy densities are per transverse area and per longitudinal rapidity, i.e., one would divide by To to get units 
per fm 3 . The three parameters are / wn , the saturation cross section cr sat , and the normalization (dE/dy) pp . When 
two identical columns of nuclei collide, T A = Tg, which leads to e wn = e sat . The quantity a nn is not an adjustable 
parameter, it is the known inelastic nucleon- nucleon cross section of 42 mb. 

In the diffuse limit, where T^,Te — >• 0, the energy density becomes (dE / ' dy) pp TATB& nni which is known as the 
binary collision limit. If one considers two diffuse nuclei colliding randomly over a large area 5, one finds the average 
energy per area in either expression to be 



(dE/drj) 



(j nn {dE/dy) T 



dxdy 



T A (x,y) J dx'dy' T B (x' \y f ) 



ABa nn (dE/dy) p 



(5) 



The parameter (dE/dy) pp is the energy per unit ra- 
pidity of a pp collision at To- Although that number is 
measured in the asymptotic limit, it might be different 
for to = 0.8 fm/c. Thus, it is a parameter that is ad- 
justed from 0.85 to 1.2 times the energy per rapidity of 
a pp collision of [25] . 

The parameter cr sat controls the scale for changing the 
behavior of e sat from the binary collision limit where 
6 ~ T A T B to the saturated limit when c ~ T m i n . The 
change occurs for T max w l/cr sat . The parameter cr sat 
also changes the wounded nucleon scaling form from that 
of binary collisions to the saturated limit where it is pro- 
portional to Ta +Tb. 



The wounded nucleon and saturation expressions dif- 
fer when T a ^ T5. For the case where cr sat T a ^> 1 and 

Vs&tTb <C 1, 

lim e wn = {dE/dy) ^ T a /2, (6) 

(JsatTa^l ^"sat 

cr sat T b <l 

lim e sat = {dE/dy)ppann 2T b . 

cr sa tT a >l (7 sat 
cr sat T b <l 

For a single nucleon, cr sat T5 <C 1, colliding onto a thick 
target, cr sat T a ^> 1, the energy density in the wounded 
nucleon expression continues to scale proportional to T a . 
For example, colliding a single nucleon onto a target with 
^sat^a = 10 6 would give nearly 1000 times the multiplic- 
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ity for a collision with cr sat T a = 1000. In contrast, the 
saturation formula would give roughly the same energy 
density for both instances. It was shown in [26] that 
differences such as these significantly affect the initial 
elliptic anisotropy, and therefore significantly affect the 
measured elliptic flow. This can be understood by con- 
sidering the collision of two equal mass nuclei with an 
impact parameter in the x direction. Along the x = 
line, both the wounded nucleon and saturation expres- 
sions give the same energy density. However, if one goes 
outward so that x becomes sufficiently large that one is 
at the edge of one nucleus, while being near the center of 
the other nucleus, the wounded nucleon formula gives a 
significantly higher energy density. This gives a relatively 
lower elliptic anisotropy for the wounded-nucleon model, 
and results in lower elliptic flow for the wounded-nucleon 
form than for the saturation form. 

The fourth varied parameter describes the initial trans- 
verse flow, i.e., the collective flow at To = 0.8 tm/x. Ini- 
tial flow has been found to significantly affect femtoscopic 
source sizes [27] and elliptic flow [28] . In [29] [30] it was 
shown that one can express the transverse flow as 
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(7) 



given four conditions: (a) A traceless stress-energy ten- 
sor, (b) Lowest order in r, (c) Bjorken boost-invariance, 
(d) Anisotropy of the stress-energy tensor independent 
of x and y. The power of the parameterization is that 
in the high-energy limit one expects each of these con- 
ditions to be reasonably met. However, at finite energy 
and for higher orders in r, (7|) can only serve as a guide 
to set a scale for the initial flow and cannot be trusted to 
better than a factor of two. For that reason, the initial 
flow is parameterized as a constant Ffl ow multiplied by 
the amount given in Q for Toi/Too- The fraction Ffl ow 
was varied from 0.25 to 1.25. 

For this first study, the initial energy density profiles 
are calculated from the average areal densities of the in- 
coming nuclei, and are smooth, as if many events from 
the same impact parameter were averaged together. This 
is known to be fairly unrealistic, and the shortcoming will 
be addressed in the future. 



B. Hydrodynamic Module 

Viscous hydrodynamics in an environment where there 
are no net conserved charges is based on local energy mo- 
mentum conservation plus two assumptions. First, it is 
assumed that in the rest frame of the stress-energy tensor 
the effective pressure equals the equilibrated pressure, 



(8) 



i.e. the bulk viscosity is assumed to be zero. Second, it is 
assumed that the remainder of the stress-energy tensor is 



sufficiently close its Navier Stokes value that its evolution 
can be described with Israel- Stewart equations of motion, 
which in the frame of the fluid becomes 



7Tj, 
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(9) 



-^V- v 



The Israel- Stewart relaxation time was set to, tis = 
3r]/sT, which is the correct expression for an equilibrated 
gas of massless particles. Once these conditions are met, 
applying the local conservation of energy and momen- 
tum, 



0, 



(10) 



determines the evolution of the stress-energy tensor. 

At high energy density the first assumption, that 
Tu = 3P(e), can be met even if the system is far from 
chemical or kinetic equilibrium. For a gas of weakly inter- 
acting massless particles, or even for a region dominated 
by weakly interacting classical fields, the condition is met 
regardless of the configuration of either the particles or 
the fields. Once the fireball cools down near the transi- 
tion region, and conformal invariance is lost, this assump- 
tion becomes questionable. The second assumption may 
be poorly met during the first 1-2 fm/c. However, the 
impact of changing the anisotropy of the stress-energy 
tensor at early times tends to be rather small [29 . 

The hydrodynamic module used here is built on an as- 
sumption of longitudinal boost-invariance which allows 
the calculations to become effectively two-dimensional 
before solving Israel- Stewart equations of motion. This 
approach has been applied by numerous research groups 
[3TH34] [62], The reduction of the dimensionality is justi- 
fied to better than the five percent level [35]. The equa- 
tion of state, P(e), comes from lattice calculations of 
Wuppertal-Budapest group [15] for temperatures above 
the hadronization temperature, and use a hadron-gas 
equation of state at lower temperatures. The equation of 
state for temperatures just above the hadronization tem- 
perature is slightly modified from the lattice values to 
match the hadron gas value at the hadronization thresh- 
old. 

For temperatures above 170 MeV/c, the viscosity to 
entropy density ratio was described with two parameters, 



V 

s 



tin 



(11) 



where T c is assumed to be 170 MeV. The first parameter, 
r]/s | t c , describes the viscosity just above the hadroniza- 
tion threshold, while the second parameter, a, describes 
the temperature dependence. 

The hydrodynamic/cascade interface temperature was 
set at a temperature of 170 MeV. Calculations were 
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also performed for a hadronization temperature of 155 
MeV, but those calculations consistently over-predicted 
the flow, or equivalently, under-predicted the number of 
hadrons for a given amount of transverse energy. It is 
the author's intention to perform a detailed study of the 
sensitivity to the equation of state and the details of 
hadronization in a separate paper. A summary of model 
parameters is provided in Table. [l[ 



C. Hadronic Cascade 

The hydrodynamic module was run until all elements 
cooled below 170 MeV. During the hydrodynamic evo- 
lution, the properties of the 170 MeV hypersurface were 
recorded. This included the boundaries, the flow velocity 
and the anisotropy of the stress-energy tensor. Hadrons 
were generated with a Monte Carlo procedure ensuring 
that all elements of the stress-energy tensor were contin- 
uous across the hyper-surface. The method [36] assumes 
that all species have a single time relaxation scale in- 
dependent of their momentum. Other approaches have 
considered the effect of adding a momentum or species 
dependence to the relaxation time [32 , but because this 
study considers only particles with low to moderate p t , 
and because the particles interact a few more times in 
the cascade module, the details of the algorithmic choice 
are not expected to matter, as long as the stress-energy 
tensor is continuous across the boundary. 

A list of particles as produced in the interface was then 
fed into the cascade on an event-by-event basis. For these 
studies, 4000 cascade events were produced for each im- 
pact parameter. The cascade code was inspired by the 
physics of the hadronic module of URQMD [37 , but was 
significantly rewritten avto improve speed, and is labeled 
B3D [58] . Hadrons were assumed to collide through res- 



onances with Breit-Wigner forms, plus a simple s— wave 
elastic cross section of 10 mb. The s— wave cross sec- 
tion was chosen independent of momentum and parti- 
cle species. The resonances from the particle data book 
[39] with masses less than 2.2 GeV/c 2 were all included. 
On average, particles collided less than twice after be- 
ing generated from the hydrodynamic interface. Pions 
had fewer collisions on average, while protons had more. 
The collisions in the cascade mainly affected the spec- 
tra and V2 of protons. There are numerous ways to im- 
prove the cascade, such as more realistic cross sections, 
consistent time delays in scattering processes, mean-field 
effects, and Bose effects for pions. However, given the 
rather modest impact of the cascade at high energy, it is 
not expected that the observables would change signifi- 
cantly. 

The B3D code runs approximately two orders of mag- 
nitude faster than URQMD for the calculations used 
here. This is mainly due to two improvements: bet- 
ter handling of the linked lists used to track collisions, 
and adding cyclic boundary conditions so that boost- 
invariance could be efficiently incorporated. The major- 
ity of the numerical expense of the calculations came from 
the cascade, and improving the speed allowed a greater 
number of points in parameter space to be explored. 

The cascade ran until all collisions ceased. For 
each outgoing particle, the momentum, particle id and 
the space-time coordinates of the last interaction were 
recorded. Since the reaction plane is known, it is straight- 
forward to calculate the azimuthal anisotropy factor v 2 = 
(cos 2(j)) . Spectra are efficiently calculated given that the 
cyclic boundary conditions make it possible to use all the 
particles when calculating the spectra at zero rapidity. 

D. Final-State Interactions 



Two-particle correlations at small relative momentum provide femtoscopic information about the phase space dis- 
tributions. This information is expressed through the Koonin formula [40] |4T|. 

C(K = ( Pl + p 2 )/2, k = ( Pl - p 2 )/2) = J d 3 r 5(K, r) |0 k (r) | 2 , (12) 

_ / d 3 M 3 r 2 /(K, n)/(K, n)5(r - (n - r 2 )) 



S(K,r) 



f<P ri <Pr 2 /(K,n)/(K,n) 



Here, </> q (r) is the outgoing two-particle wavefunction, 
/(p, r) is the phase space density in the asymptotic state, 
and *S(K, r) describes the chance that two particles with 
the same asymptotic momentum K would be separated 
by r should they not interact. Correlations provide 
the means to determine the coordinate-space informa- 
tion of 5(K, r) from the measured correlations, C(K, k). 
Through a fitting procedure, one can infer source radii 
which fit the shape of *S(K, r) with Gaussian radii, i.e. 



5(K,r) ~ expi-a^L _ y y 2 R^ e - z 2 /2B? ong }, 
where the "outward" direction is transverse to the beam 
and parallel to K, the "longitudinal" direction is along 
the beam axis and the "sideward" direction is perpendic- 
ular to the other two. 

The source radii are typically extracted by experimen- 
tal collaborations through fitting their measured correla- 
tions to expectations from Gaussian sources. Descrip- 
tion of such analyses can be found in [41] . For the 



6 



parameter 


description 


range 


(dE/dy) pp 


The initial energy per rapidity in the diffuse limit compared to measured value in pp collision 


0.85-1.2 


C"sat 


This controls how saturation sets in as function of areal density of the target or projectile. In 
the wounded nucleon model it is assumed to be the free nucleon-nucleon cross section of 42 mb 


30 mb-50 mb 


fwn 


Determines the relative weight of the wounded-nucleon and saturation formulas for the initial 

energy density described in (|3| |4| 


0-1 


•fflow 


Describes the strength of the initial flow as a fraction of the amount described in ( 7 ) 


0.25-1.25 


v/s\t c 


Viscosity to entropy ratio for T = 170 MeV 


- 0.5 


a 


Temperature dependence of rj/s for temperatures above 170 MeV/c, i.e., 
ri/s = r 1 /s\ Tc +aln(T/T c ) 


- 5 



TABLE I. Summary of model parameters. Six model parameters were varied. The first four describe the initial state being 
fed into the hydrodynamic module, and the last two describe the viscosity and its energy dependence. 



model calculations correlation functions were calculated 
by first sampling S'(K, r) then combining pairs of pions 
with similar momentum. Pions were divided into bins 
of 20 MeV/c width in transverse momentum and in 15° 
bins in azimuthal angle, before pairing. Utilizing boost 
invariance, all the pions could be longitudinally boosted 
to a frame where the rapidity was zero. The space-time 
points at which particle's had their last interaction had 
been recorded along with their asymptotic momentum 
during the running of the B3D module. This allowed a 
list of r = i*i — T2 to be constructed for each momen- 
tum bin. Correlation functions for each momentum bin 
were calculated by assuming a simplified wave function, 
|0 q (r)| 2 = l + cos(2k-r). Gaussian source radii were then 
found by searching for radii that best reproduce the cor- 
relation functions calculated by the model. Thus, rather 
than matching experimental and theoretical correlation 
functions, Gaussian radii were compared. The calcula- 
tion of correlation functions and fitting was performed 
with the code base in CorAL [42] . 



III. REDUCTION OF EXPERIMENTAL DATA 
FOR STATISTICAL ANALYSIS 

The heavy ion data sets from RHIC and from the 
Pb+Pb experiments at the LHC represent some of the 
largest scientific data sets in existence. A principal moti- 
vation of this work is to develop a statistical analysis that 
can be extended to large heterogenous data sets. This 
would include data taken at multiple beam energies, with 
different target-projectile combinations and with differ- 
ent detectors. The recent beam-energy scan at RHIC and 
the inauguration of the LHC have increased the available 
data by more than an order of magnitude as compared to 
the Au+Au collisions at 100A GeV beams measured at 
RHIC. Additionally, analyzed measurements of Cu+Cu, 
Cu+Au and U+U from RHIC will soon be available. The 
data set from the one beam energy contains petabytes of 
information. For this first study, we confine our analysis 
to this one data set, Au+Au at 100 A GeV + 100 A GeV. 
We further confine the analysis to a subset of soft physics 
observables: spectra, elliptic anisotropy, and femtoscopic 
correlations. Only mid-rapidity observables were consid- 



ered. These are the observables most connected to the 
bulk dynamics and to the bulk properties of matter, and 
are often referred to as "soft physics" . Several classes of 
observables are being ignored, e.g., jet quenching, long- 
range fluctuations and correlations, dilepton and direct- 
photon measurements, and heavy flavor. These observ- 
ables are often labeled "rare probes" and their interpre- 
tation largely factorizes out of the analysis of the soft ob- 
servables being considered here. For instance, although 
jet quenching depends on the energy density and bulk 
properties of the quark gluon plasma, the soft physics 
observables being considered here are not significantly af- 
fected by the mechanism for jet production. Further, the 
theory and phenomenology governing these other classes 
of observables often carry large uncertainties, not only 
in additional unknown parameters, but also in that they 
carry questions concerning the choice of approach. Given 
the way that the physics from these other classes of anal- 
yses factorize from the soft physics, and the lack of the- 
oretical consensus, the prudent course of action seems 
to be to determine the bulk dynamics of the system us- 
ing the soft physics observables. Once the evolution of 
the system is determined, with quantified uncertainties, 
one would have a validated basis from which to calculate 
other classes of observables, such as rare probes. 

Within the set of soft-physics observables, this first 
analysis is restricted to a subset of the overall data. 
For spectra, we consider only pions, kaons and protons. 
It would be straight-forward to consider strange baryon 
spectra, but due to large systematic and statistical errors, 
they are unlikely to greatly affect the answer at the cur- 
rent time. Additionally, because theoretical treatments 
away from mid-rapidity remain in an immature stage, our 
analysis concerns only mid-rapidity observables. For an- 
gular anisotropics, we consider only V2 and ignore higher 
order anisotropics for n > 2, 

v n = (cos(n0)), (13) 

where <j) is the angle of a particle relative to the reac- 
tion plane. Recent analyses of v n> 2 suggest that the 
observables may even be more sensitive to the viscosity 
than V2 [43] [44] . However, theoretical questions remain 
about how to instantiate the event-by-event fluctuations 
which drive these higher-order harmonics. This analysis 



7 



only considers V2 for pions. Although v 2 is measured for 
kaons and protons, in order to compare to data, theoreti- 
cal treatments would have to run for tens of thousands of 
events for each impact parameter to get sufficient statis- 
tics for kaons and protons. This analysis used 4000 events 
per impact parameter. Finally, the femtoscopic analysis 
is confined to same-sign pions. Source sizes extracted 
from other analyses carry significantly more uncertainty. 
RHIC data is recorded according to centrality bins, e.g., 
top 5%, top 10%, 0-10% — Bins are typically assigned 
according to some measure of overall multiplicity. For in- 
stance, the 20-30% bin corresponds to those events with 
multiplicities that are lower than the top 20% of events 
and higher than the lower 30% of events. The choice of 
bins varies between observables and between collabora- 
tions. Since the hydrodynamic treatment is questionable 
at low centrality, we decided to neglect the more periph- 
eral collisions, i.e. those with centrality greater than 30%. 
Further, analysis was confined to two bins, 0-5% and 20- 
30%, due to the expectation that if those two bins were 
matched, any intermediate bin would also be matched. 
This reduced the numerical cost of performing the simu- 
lations. 

It is our hope to extend future analyses to include 
more data. This would include data from the RHIC 
beam-energy scan, from the LHC, and from the Cu+Cu, 
Cu+Au and U+U at RHIC. Data from the LHC is 
straight-forward to incorporate because the same theo- 
retical models can be used once one has added energy- 
dependence to the initial-state parameterization. The 
Cu+Au and U+U data require significantly rethinking 
the parameterization of the initial state, especially for 
Uranium due to large nuclear deformation. Extending 
the analysis to include data from the beam-energy scan 
would require significant changes to the model used here. 
At lower energies, one can no longer assume Bjorken 
boost-invariance and can no longer ignore the baryon 
excess. Although our present hydrodynamic code can 
work in three dimensions, significant theoretical work 
is required to develop a parameterization for the three- 
dimensional initial state at arbitrary beam energies. 



A. Initial Distillation of Observables 

Experimental collaborations spend a lot of effort re- 
ducing the huge RHIC data set to a finite number of 
published plots representing a useful summary of some 
observable. For example, the PHENIX and STAR col- 
laborations have produced plots of proton spectra for 
several centrality classes. Each plot might have a few 
dozen points. It is infeasible for our analysis to consider 
reproducing each of these data points. Instead, each ob- 
servable was reduced to a few represenative quantities. 
For a given species and centrality, spectra were reduced 
to two numbers. The first number is the yield, or inte- 
grated spectra, within a finite p t (transverse momentum) 
range. The ranges were set so that they ignored the high 



p t tail, which is strongly affected by jets and is outside 
the scope of the model. The second number would be 
the mean p t within that range. The choice to use mean 
p t was justified by doing a principal component analysis 
(PCA, see below) on the data points within a spectra di- 
vided by the yield. This ensured that only the shape of 
the spectra was coming into play. 

To some extent such an analysis is inherently arbitrary, 
one needs to understand both theoretical and experimen- 
tal systematic errors. By "theoretical systematic errors" , 
we refer to the fact that even if all the model parameters 
were correct, shortcomings of the model would preclude 
one from making comparing to data beyond a certain ac- 
curacy. Both classes of systematic error could have off- 
diagonal elements if one were to construct aniVxiV error 
matrix, where N would be the number of data points in 
a specific spectra. The first principal component is sen- 
sitive to exactly how such a matrix would be generated. 
After trying several competing choices for systematic er- 
ror, it was clear that the shape of the spectra was well 
represented by a single component. Further, that com- 
ponent was similar to the mean p t . In essence, any two 
model runs that yield two spectra with the same p t also 
yield indistinguishable spectra up to uncertainties. Thus, 
the spectra can be reasonably well summarized by the 
yield and mean p t . 

For the elliptic flow, the experimental information con- 
sisted of plots of V2 as a function of p t . A PCA analysis 
showed that the pt-weighted value for V2 effectively cap- 
tured all the information within the set of model runs. 
Femtoscopic information was contained by plots from the 
STAR Collaboration of Gaussian radii (i? Q ut? ^side and 
^long) as a function of the transverse momentum. Simply 
averaging each radius over the several p t bins was found 
to effectively encapsulate nearly all the variation of the 
femtoscopic radii throughout the model runs. 

In this manner the various experimental results were 
thus reduced to those listed in Table. Till Each observable 
was also assigned an uncertainty. This uncertainty rep- 
resented the accuracy within which a comparison of the 
theoretically determined value from a model run could be 
meaningfully compared to the corresponding experimen- 
tal measurement. Of all the observables in Table. [TT] only 
V2 has significant statistical error. Additionally, the 
observable is known to be significantly affected by known 
shortcomings in the model, such as the lack of event-by- 
even fluctuations. Thus, V2 is assigned a larger percent- 
age error than other observables. 



B. PCA Analysis of reduced observables 

One could create model emulators for each of the ob- 
servables listed in Table, [ill However, one can further dis- 
till the data to a handful of principal components repre- 
senting their most informative linear combinations. This 
serves to further reduce the complexity of the surrogate 
model. Let y exPi i and <ii be data points and uncertainties 
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observable 


pt weighting 


centrality 


collaboration 


uncertainty 


^2,7r + 7r- 


ave. over 


11 p t bins from 160 MeV/c to 1 GeV/c 


20-30% 


STAR* g5] 


12% 


Rout 


ave. 


over 4 p t bins from 150-500 MeV/c 


0-5% 


STAR g6] 


6% 


Rside 


ave. 


over 4 bins from 150-500 MeV/c 


0-5% 


STAR [46] 


6% 


-Rlong 


ave. 


over 4 pt bins from 150-500 MeV/c 


0-5% 


STAR [46] 


6% 


-Rout 


ave. 


over 4 bins from 150-500 MeV/c 


20-30% 


STAR [46] 


6% 


-Rside 


ave. 


over 4 pt bins from 150-500 MeV/c 


20-30% 


STAR [46] 


6% 


Rlong 


ave. 


over 4 p t bins from 150-500 MeV/c 


20-30% 


STAR g6] 


6% 






200 MeV/c < pt < 1.0 GeV/c 


0-5% 


PHENIX [47] 


3% 


(Pt) K + K- 




400 MeV/c < p t < 1.3 GeV/c 


0-5% 


PHENIX [47] 


3% 


(Pt)pp 




600 MeV/c < pt < 1.6 GeV/c 


0-5% 


PHENIX [47] 


3% 






200 MeV/c < pt < 1.0 GeV/c 


20-30% 


PHENIX [47] 


3% 


(Pt) K +K- 




400 MeV/c < pt < 1.3 GeV/c 


20-30% 


PHENIX [47] 


3% 


(Pt)pp 




600 MeV/c < pt < 1.6 GeV/c 


20-30% 


PHENIX g7] 


3% 


7v + tv~ yield 




200 MeV/c < pt < 1.0 GeV/c 


0-5% 


PHENIX [47] 


6% 


7r + 7r _ yield 




200 MeV/c < pt < 1.0 GeV/c 


20-30% 


PHENIX [47] 


6% 



TABLE II. Observables used to compare models to data. *To account for non-flow correlations, the value of V2 was reduced 
by 10% from the value reported in [45] . 



for the i = l through N data points listed in Table. 
One then considers the corresponding quantities from the 
model run m, y m ^ where m runs from 1 to the number 
of full model runs M. A useful first step is to scale the 
quantities by their net uncertainty, 

-(Vi) 



Ve 



Ve 



Vn 



(14) 



(Vi) 



1 M 

m=l 

The net uncertainties, G{ are operationally defined as the 
uncertainty involved in comparing a model value to an 
experimental measurement. The measurements consid- 
ered in this paper are limited by systematic rather than 
aleatoric errors, and we assume that errors are described 
by a normal distribution, 



£(x) ~ exp < 



O ex p) _ 7/ ( mod ) 



to) 



2 \ 



2a 



(15) 



where ?/ exp ) and ?/ mod ) are the experimentally measured 
and model values respectively. Even if the model param- 
eters are exact, the models also have limited accuracy 
due to shortcomings in the physics. Thus, the net un- 
certainty encapsulates both theoretical and experimental 
uncertainties, i.e., they can be considered to describe the 
inability of the model not only to describe the physics of 
the collision, but to also account for the inadequacy of 
the model to describe uncertainties in the experimental 
measurement and analysis. 

One then studies the sample covariance of the model 
values amongst the M model runs. 



1 M 



(16) 



The N eigenvalues of S are A^, and the normalized eigen- 
vectors are €ij. One can then consider new variables, z m ^ 
which are linear projections of the original y m ^ along the 
various directions defined by the eigenvectors, 



3 



(17) 



With this procedure, the model values, y m i, are rotated 



into a basis where the values 
variance over the model runs, 



have a diagonalized 



M 

-Y 

If ^ 



M 



(18) 



rn=l 



The values z m ^ are known as principal components. 
Since the values y were scaled by the uncertainties, the 
components yi have uncertainties of unity, and after rota- 
tion the values Zi also have uncertainties of unity. Since 
the variance of z within the model runs is diagonal, one 
can state that those components for which <C 1 can 
be ignored because they do not assist in discriminating 
parameters. Further, the discriminating power is often 
dominated by the first few principle components, i.e., 
those with the largest A^. 

To further justify our selection of principle components 
we show a plot of the normalized cumulative variance 
explained by the largest r components in Figure, [l] i.e 



F(r) 



(19) 



where we have sorted the eigenvalues into descending or- 
der. Examination of this figure clearly shows that the 
first four principle components are sufficient to explain 
almost all of the sampled variance. 

Once the principal components, Zi, have been deter- 
mined, one can invert the transformations to find yi in 
terms of the Z{. The components which do not contribute 
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5 10 

No. of Principal Comp.s, r 



FIG. 1. The variance resolving power R(r) of the princi- 
ple components, only the first few components are needed to 
explain almost all of the observed variance. 



strongly to the total variance can be set to zero and the 
resulting ^s will not be appreciably affected. In this par- 
ticular case these are the components with A< 1. Thus, 
the statistical analysis need only emulate those compo- 
nents with Xi > 1. 

Given the 15 observables outlined in Table. |TTJ one 
could construct an emulator for each observable. How- 
ever, a PCA analysis of the 15 intermediate observables 
shows that not more than four of the principal compo- 
nents vary appreciably throughout the model runs. For 
these four components, the corresponding fluctuations, 
(Szf) were of order unity or greater, while the remain- 
ing components fluctuated much less than unity. Thus, 
instead of tuning 15 emulators, only six principal compo- 
nents were considered (even though only four were truly 
needed). It is instructive to list the values of A^ and the 
decomposition of the main components. This is shown 
in Table. |III| The eigenvalues A^ represent the resolving 
power of the various principal components. 

To gain an understanding of the degree to which the 
parameter space is constrained by each measurement of 
Zi, one can consider the case where the prior distribution 
of parameters, x ai is distributed with unit variance ac- 
cording to a Gaussian distribution, (x a xp) = 5 a p. We 
further assume that each principal component Zj varies 
linearly with x. In that case, the gradient of each princi- 
pal component, 



dz ■ 



(20) 



form a set of orthogonal vectors because the covariance 
(ziZj) = XiSij is diagonal, 

(ziZj) = (Vzi) a (Vzj)p(x a xp) = Vzi ■ Vzj = XSij (21) 

Thus, if each component of z depends linearly on x, each 
principal component constrains a separate direction in 
parameter space. One can then understand the resolving 



power by considering the simple case with one principal 
component and one parameter. Given a measurement 
^(exp) an d assuming that the prior has unit variance and 
that z depends linearly with x, 



P(x) ~e ■ 1 ' 2/2 exp <j - (rnx - z (exp) 



) 2 /2}, (22) 



where dz/dx is the slope m, and from (21 ) m 2 = A. Com- 
pleting the squares in the argument of the exponential, 



P(x) 



-(X+l)(x-n 2 )/2 



(23) 



where \i is the prior mean for x. This shows that if the 
response is purely linear that each principal component 
reduces the width of the posterior relative to the prior by 
a factor 1 / y/l + A$. 

The first principal component in Table. [TTT| carries the 
bulk of the resolving power. Since Ai = 18.36, the linear 
considerations above suggest that a measurement of the 
first principal component should constrain the original 
parameter space by a factor of roug hly 1/Vl9l. The 
second and third principal components also significantly 
narrow the parameter space. All together, one expects 
these measures to constrain the fit to on the order of 5% 
of the original six-dimensional parameter space at the 
"one-sigma" level. This estimate of the resolving power 
is based on an assumption that z varies linearly with x, 
but nonetheless provides a useful, although crude, ex- 
pectation for how the our analysis might ultimately con- 
strain the parameter space. 

The first and second components dominantly consist 
of measures of the multiplicity and of the V2 observable. 
This is not surprising. It shows that the most important 
aspect of fitting data is to fit the multiplicity and elliptic 
flow. The third component has a large mixture of (p t ) 
and inter ferometric observables. Thus, before performing 
the parameter space exploration, one expects that those 
parameters driving the multiplicity and elliptic flow will 
be the most significantly constrained. 



IV. THEORY OF MODEL EMULATORS 

Developing an understanding of a six-dimensional pa- 
rameter space requires many millions of samplings of the 
model output. Running a complex code for each sam- 
pled point in parameter space is impractical. An alterna- 
tive strategy has been to develop model emulators, codes 
which effectively interpolate from an initial sampling of 
runs through the space . 

We construct a Gaussian Process emulator [48H5T] , 
which acts as a statistical model of our computer model. 
An emulator is constructed by conditioning a prior Gaus- 
sian Process on a finite set of observations of model out- 
put, taken at points dispersed throughout the param- 
eter space. Once the emulator is trained it can rapidly 
give predictions for both model outputs and an attendant 
measure of uncertainty about these outputs at any point 
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observable \ A* 


18.36 


7.87 


0.93 


0.21 


0.04 


0.012 


cent0to5_PHENIX_spectraPION_YIELD 


0.43202 





52170 


0.21636 


0.56290 


0.06883 


0.35417 


cent0to5_PHENIX_spectraPION_MEANPT 


0.10117 





02647 


0.37032 


-0.08869 


0.09235 


-0.24640 


cent0to5_PHENIX_spectraKAON_MEANPT 


0.10770 





03291 


0.37755 


-0.07459 


0.06328 


-0.26766 


cent0to5_PHENIX_spectraPPBAR_MEANPT 


0.04925 





02192 


0.16751 


-0.02131 


-0.05466 


-0.19039 


cent0to5_STAR_ROUT_PION 


-0.01942 





06908 


-0.31734 


-0.02968 


0.72626 


0.12886 


cent0to5_STAR_RSIDE_PION 


0.09148 





09321 


0.07972 


0.05565 


0.11943 


-0.07137 


cent0to5_STAR_RLONG_PION 


0.08413 





09520 


-0.13599 


0.37546 


0.08521 


-0.50343 


cent20to30_PHENIX_spectraPION_YIELD 


0.43743 





49869 


-0.32721 


-0.56043 


-0.26805 


-0.01286 


cent20to30_PHENIX_spectraPION_MEANPT 


0.07549 





03028 


0.33981 


-0.23142 


0.28313 


-0.06472 


cent20to30_PHENIX_spectraKAON_MEANPT 


0.08266 





03721 


0.34043 


-0.23645 


0.27941 


-0.06785 


cent20to30_PHENIX_spectraPPBAR_MEANPT 


0.03791 





02697 


0.14297 


-0.11517 


0.03747 


-0.09339 


cent20to30_STAR_V2_PION_PTWEIGHT 


-0.74299 





65843 


0.08846 


-0.03607 


-0.01531 


-0.06192 


cent20to30_STAR_ROUT_PION 


0.02955 





03296 


-0.30420 


-0.06375 


0.43249 


-0.09820 


cent20to30_STAR_RSIDE_PION 


0.08368 





09367 


-0.01379 


-0.21381 


0.08021 


-0.06598 


cent20to30_STAR_RLONG_PION 


0.08974 





08905 


-0.24592 


0.19088 


-0.07458 


-0.62873 



TABLE III. The first six principle components. Since the variables were initially scaled by their uncertainties, the eigenvalues, 
Ai, describe the resolving power of the components. Only the first ~ 4 components are significant, i.e., A > 1. The table also 
provides the decomposition of the principal components in terms of the 15 observables. 



in the parameter space. This is a probability distribution 
for the model output at all points in parameter space 
and is by far the useful feature of Gaussian Process em- 
ulators. The most common interpolation schemes, such 
as interpolating polynomials, produce an estimate of the 
model output at a given location in the parameter space 
with no indication as to the extent that this value should 
be trusted. Furthermore, numerical implementations of 
Gaussian Process emulators are computationally efficient 
(producing output in microseconds rather than minutes), 
making it feasible to predict vast numbers of model out- 
puts in a short period of time. This ability opens new 
doors for the analysis of computer codes which would 
otherwise require unacceptable amounts of time [52j [53] . 

We construct an emulator for a model by conditioning 
a Gaussian Process prior (see Figure. [2| on the train- 
ing data [54H56] . A Gaussian Process is a stochastic 
process with the property that any finite set of sam- 
ples drawn at different points of its domain will have a 
multivariate-normal (MVN) distribution. Samples drawn 
from a stochastic process will be functions indexed by a 
continuous variable (such as a position or time) as op- 
posed to a collection of values as generated by, e.g., a 
normally-distributed random variable. A Gaussian Pro- 
cess is completely specified in terms of a mean and co- 
variance, both of which can be functions of the index- 
ing variable x. The covariance, c(xi,X2), might be any 
positive- definite function of x\ and £2- An example of 
unconditioned draws is shown in the left panel of Fig- 
ure. [2] for the case where the covariance depends only on 
x\ — X2 and is a power-exponential covariance function 
with unit length. The draws are smooth functions over 
the domain space, and if enough samples are drawn from 
the process the average of the resulting curves at each 
point would converge to zero. 

A predictive distribution for the value of a computer 
model at new points in the design space can be obtained 
by conditioning this process on a set of training points 
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FIG. 2. Left panel: Unconditioned draws from a Gaussian 
Process GP(0, 1) with a mean of zero and constant unit vari- 
ance. Right panel: draws from the same process after condi- 
tioning on 7 training points (black circles) . The gray band in 
both panels is a pointwise 95% confidence interval. Note how 
the uncertainty in the right panel grows when away from the 
training points. 



obtained from running the model. Conditioning forces 
samples drawn from the process to always pass through 
the training points. The resulting curves interpolate the 
training data, as shown in the right hand panel of Fig- 
ure. [2] Repeated draws from the conditioned posterior 
distribution would on average follow the underlying curve 
with some variation, shown by the gray confidence re- 
gions. These confidence bubbles grow away from the 
training points, where the interpolation is least certain, 
and contract to zero at the training points where the 
interpolation is absolutely certain. The posterior distri- 
bution can be evaluated to give a mean and variance at 
any point in the parameter space. We may interpret the 
mean of the emulator as the predicted value at a point, 
the variance at this point gives an indication of how close 
the mean can be expected to be to the true value of the 
model. 
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To construct an emulator we need to fully specify our 
Gaussian Process (GP) by choosing a prior mean and a 
form for the covariance function. The model parameter 
space is taken to be p-dimensional. We model the prior 
mean by linear regression with some basis of functions 



h(x). In this analysis we use the trivial basis h(x) = {1}. 
We specify a power exponential form for the covariance 
function with power a ~ 2 to ensure smoothness of the 
GP draws (a has to be in [1,2] to ensure positive defi- 
niteness), 



c(xi,Xj) = o exp [-^2 



k=l 



Qk 



SijO N , ae[l,2]. 



(24) 



Here, #o is the overall variance, the k set characteristic 
length scales in each dimension in the parameter space 
and On is a small term, usually called a nugget, added to 
ensure numerical convergence or to model some measure- 
ment error in the code output. The shape of the covari- 
ance function sets how the correlations between pairs of 
outputs vary as the distance between them in the param- 
eter space increases. The scales in the covariance function 
k are estimated from the data using maximum likelihood 
methods [56], in Figure. [3] we demonstrate their influence 
on an artificial data set. The linear regression model 
handles large scale trends of the model under study, and 
the Gaussian Process covariance structure captures the 
residual variations. 

Given a set of n design points V = {a?i, . . . , x n } in 
a p-dimensional parameter space, and a set of n train- 
ing values representing the model output at the design 
locations Y = {yi, . . . , y n } , the posterior distribution 
defining our emulator is 

P(x, 0) ~ GP (m(x, 0), E(aj, (9)) , (25) 

for conditional mean rh and covariance E. 

rh{x) = h(x) T (3 + fe T (a;)C- 1 (Y P - H/3), 
E(a?i, Xj) = c(xi, Xj) — k T (xi)C~ l k(x j) + T(xi,Xj), 
Cij = c(xi,Xj) (26) 

r(x hXj ) = (h(xi) T - k 1 \x i )C- 1 U) T (l^C -1 !!)" 1 

k(x) T = (c(x!,x), . . . , c(x n , x)) . (27) 

Where rh(x) is the posterior mean at cc, fl(xi,Xj) is 
the posterior covariance between points X{ and £Cj, C 
is the n x n covariance matrix of the design D, /3 are the 
maximum-likelihood estimated regression coefficients, h 
the basis of regression functions and H the matrix of 
these functions evaluated at the training points. 

The elements of the vector k(x) are the covariance of 
an output at x and each element of the training set. It 
is through this vector k(x) that the emulator "feels out" 
how correlated an output at x is with the training set 
and thus how similar the emulated mean should be to 
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FIG. 3. Demonstration of emulator behavior as a function 
of correlation length, 0\. In all panels, the solid blue line 
shows the mean of the emulator and the solid gray region is a 
95% confidence interval around this region. Left panel: fitting 
with a value of 6\ that is too small (under-smoothing) . Right 
panel: shows over-smoothing by using a value of Q\ that is too 
large. Central panel: smoothing with a value of 0± = 0.68 that 
was obtained by a maximum likelihood estimation method. 



the training values at those points. Note that the quan- 
tities defined in (26) depend implicitly upon the choice 
of correlation length scales = {0o,0 k ,0/v} which deter- 
mine the shape of the covariance function. 

The expression for the mean in ( |26| ) can be decom- 
posed into a contribution from the prior, the linear re- 
gression model h(x) J p plus a contribution applied to 
the residuals determined by the covariance structure 
k T (x)C~ 1 (Y — H/3). Similarly the covariance can be 
decomposed into a contribution from the prior, the co- 
variance function c(xi,Xj) plus corrections arising from 
the prior covariance structure and the covariance of the 
new location x through k(x). These terms weight the 
points Xi, Xj more highly the closer they are to the train- 
ing points through k. The Y term gives the corrections 
to the covariance arising from the regression model. 

In our study, we run the full code at N = 729 (cho- 
sen because 3 6 = 729) points from the parameter space. 
A Latin Hyper Cube (LHC) design is used to generate 
the training locations in the parameter space. This is 
an efficient design for space-filling in high dimensional 
parameter spaces. This is an efficient sparse design for 
high dimensional parameter spaces that is "space-filling" 
in the sense that all its lower-dimensional projections are 
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distributed as evenly as possible 

The output from the model code is multivariate. Al- 
though fully multivariate emulator formulations do exist 
they are challenging to implement. Instead we follow the 
now somewhat standard procedure of creating emulators 
for some decomposition of the code output see eg [52| [6Q]. 
In this case we apply a principle components decompo- 
sition to the model output and build emulators for each 
significant component as detailed above. 



V. TESTING THE EMULATOR 

The goal of this section is to investigate the reliability 
and accuracy of the Gaussian-process emulator described 
in the previous section. The emulators were tuned to full 
model calculations which were performed for 729 points 
in the six-dimensional parameter space. The points were 
chosen randomly according to a maximin Latin hyper- 
cube design [57] . The model output for each run was an- 
alyzed to determine the 15 summary observables, y^ de- 
scribed in Section. HIIl The 15 numbers were transformed 
to variables which have zero mean and uncertainties 
of unity as described in (14)) From these, six principal 



components Z{ were then constructed as the linear com- 
binations of the y with largest variance throughout the 
729 sampling runs. The emulator's task was to faith- 
fully reproduce these six "observables" as a function of 
the six-dimensional parameter space x along with pro- 
viding a reasonable measure of uncertainty. To test the 
emulator, an additional 32 runs were performed using 32 
random points in parameter space. A successful interpo- 
lation should reproduce the principal components from 
each of the 32 test runs to a given tolerence. 

Tuning the Gaussian-process emulator involved esti- 
mating the hyper-parameters described in (24). The first 
attempt at finding optimized hyper-parameters used the 
same methods of [13j|56]. However, that approach was 
not robust, and often led to inaccurate emulators, with 
the accuracy defined by how well the emulator repro- 
duced data from the 32 test runs. A more accurate result 
ensued by simply setting the hyper-radii, the O l values in 
(24), equal to half the range for each parameter x % in 
the model space. The exponent a was set to 1.5 and 
the nugget Oo was set to zero. Changing the hyper-radii 
by factors of two, or adjusting the exponent anywhere 
between 1.0 and 2.0 had little effect. For perspective, 
competing interpolating schemes were constructed, one 
based on a quadratic fit, and a second based on a linear 
fit where neighboring points were more heavily weighted 
in the fits. Each of these schemes was slightly less accu- 
rate than the Gaussian process emulator with the hyper- 
parameters chosen as described above. However, all these 
procedures performed better than the Gaussian process 
emulator using Maximum-Likelihood-Estimation (MLE) 
hyper-parameters as described in p~3j [56] . This failure to 
find good hyper-parameters may come from the numeri- 
cal challenges of the MLE optimization process given the 
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FIG. 4. (color online) The emulator error, EE, is shown for 
the 32 test runs. If the emulator error per principal compo- 
nent were of order of the experimental and model uncertain- 
ties, the values of EE would be near six. The errors above are 
significantly smaller. 



large number of training points. 

The Gaussian process emulator explicitly reproduces 
Zi(x) whenever x approaches one of the training points, 
x n ,n = 1 . . . 729. To test the emulator points had to 
be chosen away from the training points, and 32 addi- 
tional full model runs were performed at random points 
throughout the parameter space. The emulator error can 
be summarized as 



EE(x) ee £ (z\ 



(emu) 



(x) 



(mod) 



(x) 



(28) 



where z^ emu ^(x) is the conditional mean from the z'th 
emulator as given by (26) and r is the total number of 



observation principle components retained. A plot of this 
for the withheld data points is displayed in Figure. [4] for 
the 32 test runs. By construction, EE(x n ) is zero for 
the 729 training runs. The fact that the net errors were 
less than unity, even after summing over six principal 
components, shows that the emulator did an outstanding 
job of reproducing the data. Furthermore, the emulator 
error is of the order of the statistical error of the model 
(which mainly comes from the calculation of U2, which 
suggests that in this case further improving the emulator 
would not significantly improve the final result. 

The success of the Gaussian Process emulation over a 
wide range of schemes and parameters is probably due 
to the smooth and monotonic response of the model to 
parameters. At high centre of mass energies, the physical 
system is highly explosive. Within the range of param- 
eters considered the explosiveness is modified, but the 
behavior never changes qualitatively, and one expects a 
mostly linear response to the parameters. This may not 
be as true at lower energies. 

It was seen that the estimate of the errors of the emula- 



tion as defined in ( 26 ) often significantly underestimated 
the accuracy of the emulator as tested in Figure. [4] The 
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net error tended to be less than a half unit, even though 
it was summed over multiple degrees of freedom. Since 
the error associated with the accuracy of the emulator 
was so small, the emulator error was incorporated into 
the calculation of the likelihood in a simplified manner. 
The uncertainty inherent to the data and models for a 
specific principal component was unity due to the choice 
in how to scale the Zi values. By adding in the emulator 
error, the total uncertainty should be a 2 = 1 + cr 2 for 
each component. The likelihood used by the MCMC is 
then, 

^^{-i p^^f^ l (») 

For our MCMC calculations, <j e was set to 0.1 according 
to an estimate of the error per degree of freedom from 
Fig. [4] This increased the width of the posterior region 
of parameter space by only a few percent. 

VI. MCMC RESULTS 

As shown in the previous section, the emulator ac- 
curately reproduces the log-likelihood. For the MCMC 
search the Gaussian process emulator was run sampling 
many millions of points in parameter space. The trace 
provides an ergodic sample of the allowed regions in 
parameter space, i.e., the posterior distribution. The 
MCMC procedure applied here is a Metropolis algorithm. 
First, the parameter space was scaled and translated so 
that it was centered around zero, and that the flat prior 
had unit variance, i.e., it varied from — a/3 to +y/3. First, 
a random point was chosen in the six-dimensional pa- 
rameter space xi, from which one takes a random step 
to X2 = xi + Sx. The random steps £x were chosen ac- 
cording to a six-dimensional Gaussian with the step size 
in each dimension being 0.1. The likelihoods were calcu- 
lated for each point. If the likelihood £(x2) was higher 
than £(xi), the step was accepted, and if the likelihood 
was smaller, the step was accepted with the probability of 
the ratios of the two likelihoods. After the 100,000-step 
burn-in phase, the trace was stored by writing down every 
tenth point. The resulting distribution is proportional to 
the likelihood [61 and represents an ergodic sampling of 
the posterior distribution for a uniform prior. The trace 
finished when 10 6 points were written to disk. 

To evaluate the success of the emulation, 20 points 
were randomly chosen from the MCMC trace and were 
then evaluated with the full model. The observables 
used for the original analysis were then plotted for each 
of the 20 points in parameter space. Another twenty 
points were chosen randomly from the original parameter 
space, i.e. they are consistent with the flat prior distri- 
bution. Again, the observables were calculated with the 
full model for each of these points in parameter space. 
One expects the observables for each of the 20 points 
representing the MCMC trace to reasonably well match 



the experimental data, while the points chosen randomly 
from the prior distribution should lead to a wider range 
of observables, some of which should be inconsistent with 
the data. 

Comparisons of the spectra from the model runs char- 
acterizing the prior and posterior distributions are shown 
in Figure. |5j Parameters from the posterior distributions 
lead to far superior fits, to both the yields and shape 
of the spectra. From the figure, one can see that the 
spectra for heavier particles provide more discriminating 
power. This comes from the greater sensitivity to collec- 
tive flow, and emphasizes the importance of having reli- 
able measurements of proton spectra. At RHIC, STAR's 
proton spectra are warmer than those of PHENIX, and 
their estimate of the mean p t for protons is 7% higher. 
Whereas PHENIX shows the mean p t of protons staying 
steady or perhaps slightly falling with increasing central- 
ity, STAR's analysis show a rising mean p t . If the mean 
p t were indeed higher than what PHENIX reports, the 
extracted parameters should change, e.g., the initial col- 
lective flow might come out higher. 

Figure [6] shows and the femtoscopic source sizes cal- 
culated from the same representative points in parame- 
ter space for both the prior and posterior distributions. 
The MCMC is clearly successful in identifying points in 
parameter space that when run through the full model 
matched the experimental measurement. Further, given 
that the systematic uncertainty of specifying the p t av- 
eraged V2 was assumed to be 12%, the spread of V2 vs. 
p t plots appears consistent with expectations. Although 
the overall trend of the source radii were matched by 
the model, a consistent discrepancy between the data 
and model calculations using parameters from the pos- 
terior distribution is evident. At low p tl the sideward 
and longitudinal source sizes are over-predicted by more 
than 10%, which is about double the expected system- 
atic error. Over-predicting the sizes might be due to the 
assumption of boost invar iance, or perhaps the approxi- 
mations used to account for the Coulomb interaction, or 
in justifying the Koonin formula. From analyticity, one 
expects that the R ou t and R S i<±e sizes to approach one 
another as p t — >• 0. As can be seen in the upper panel 
of Figure. [6j this does not appear to be holding true in 
the data. Either the lower range of p t (200 MeV/c) is 
not sufficiently small, or an acceptance /efficiency effect 
in the detector is affecting the result. This issue should 
be resolved if femtoscopic analyses are to be applied with 
confidence at better than the 10% level. 

From the MCMC traces, the distribution of the various 
parameters and the correlations between pairs of param- 
eters are shown for the Gaussian-process emulator in Fig- 
ure. [7| The plots along the diagonal display the range of 
acceptable values for individual parameters, integrated 
over all values of the other five parameters. Although 
over 90% of the six-dimensional parameter space is elim- 
inated at the one-sigma level, the individual parameters 
are rarely constrained to less than half their initial range 
when other parameters are allowed to vary. 
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FIG. 5. (color online) Left-side panels: Pion, kaon and proton spectra from 20 model calculations where parameters are 
randomly chosen from the prior distribution. Model calculations are blue lines, and experimental data from PHENIX is shown 
as red/green circles for positive/negative charges. Results are shown for both 5% most central, and for the 20-30% centrality 
bin. Due to lack of some chemical reactions, normalizations for kaons and pions in the model were scaled by factors of 0.85 
and 0.6 respectively. Right panels: Same as left-side but with 20 model calculations where parameters were chosen randomly 
from posterior distribution as sampled by MCMC trace. 



The first four parameters ("I.C. PP NORM", "LC. 
SAT cr", "LC, W.N. FRAC" and "LC. FLOW") deter- 
mine the initial state fed into the hydro. The first pa- 
rameter "LC. PP NORM" sets the constant of propor- 
tionality between the product of the areal densities of 
the incoming nuclei, and the initial energy density fed 
into the hydro. In the limit of low aerial densities this 
should be consistent with pp collisions. Thus, the range of 
the prior distribution was quite small, and the statistical 
analysis did little to further constrain it. The parameter 
"LC. SAT a" refers to cr sat in (|3| and parameterizes the 
saturation of the energy density with multiple collisions. 
The preferred value appears rather close to the value of 42 
mb typically used in the wounded nucleon model, though 
there is a fairly wide range of accepted values. The pa- 
rameter "LC. W.N. FRAC" sets the weights between the 
wounded nucleon and the saturation parameterizations 
in ([!]). This shows a preference for the wounded nu- 
cleon prescription which gives a smaller initial anisotropy 



than the saturation parameterization. The final initial- 
condition parameterization, "I.C. FLOW" sets the initial 
transverse flow set in the hydrodynamic calculation. The 
parameter sets the initial flow as a fraction of the amount 
described by ([7]), which should be expected in the limit of 
high-energy. The MCMC trace points to a rather small 
fraction of this flow, though like all of the initial-condition 
parameters has a fairly broad range of possible values. 

The last two parameters refer to the viscosity. The 
viscosity at T = 170 MeV is referred to as "77/s" in Fig- 
ure. [7| and the temperature dependence is labelled by "T 
DEP. of 77", and refers to the parameter a in ( pT| ). Both 
are significantly constrained as a fraction of the original 
parameter space. The range of r]/s is consistent with 
similar, but less complete, searches through parameter 
space using similar models [2j [3]. In [62], the authors 
found little sensitivity to the viscosity at higher temper- 
atures, but considered a smaller variation of the viscosity 
with temperature than was considered here. 
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FIG. 6. (color online) Upper left panel: For 20 points in parameter space randomly chosen from the prior distribution, V2 
for pions is plotted as a function of p t for full model runs. Blue lines represent model calculations whereas are squares are 
experimental data. Lower left panel: Same as upper panel, except 20 points are randomly taken from posterior distribution as 
sampled by the MCMC trace. 

Right-side panels: Femtoscopic radii are shown for calculations from the prior distribution (left half) and from the posterior 
(right half) calculations. The posterior calculations well reproduce the data except for the sideward and longitudinal radii at 
low p t . 
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Figure [7] also displays cross-correlations from the 
MCMC traces. Several parameters are strongly corre- 
lated. For instance, the energy normalization "I.C. PP 
NORM" and "I.C. SAT a" are strongly correlated in 
that one can have less saturation of the cross section 
if the energy normalization is turned down. There is 
also a strong correlation between "I.C. FLOW" and "I.C. 
W.N. FRAC". One can compensate for less initial flow 
if the saturation formula is more heavily used than the 
wounded nucleon formula. Again, this is expected be- 
cause the wounded nucleon parameterization leads to less 
spatial anisotropy and a somewhat more diffuse initial 
state. 

The inferred viscosity is clearly correlated with the 
weighting between the wounded nucleon and saturation 
parameterizations, as expected from the arguments in 
[26] . The two viscous parameters are also correlated with 
one another as expected. One can compensate for a very 
low viscosity at T = 170 MeV by having the viscosity rise 
quickly with temperature. Figure [8] shows the viscosity to 
entropy ratio as a function of temperature corresponding 
to the 20 random samples from the both the prior and 
posterior distributions. Higher values of the temperature 
dependence a are increasingly unlikely for higher values 
of rj/s\ Tc . 



VII. SUMMARY AND OUTLOOK 

Two principal conclusions can be taken from this 
study. First, the data from relativistic heavy ion colli- 
sions is well suited to a multi-dimensional analysis featur- 
ing model emulators. The response of the data to model 
parameters appears sufficiently smooth to warrant sim- 
ple interpolation of a few principle components. Only 
a half dozen parameters were varied in this study, and 
only a limited number of observables were considered. 
Nonetheless, the procedure should easily scale to larger 
numbers of parameters and larger data sets. The suc- 
cess of the emulators in reproducing model output, and 
in the MCMC procedure in identifying likely regions in 
parameter space provide hope that the field can produce 
quantitative statements concerning the bulk properties 
and dynamics of the matter formed in heavy ion colli- 
sions. The second conclusion centers on the extracted 
parameters. Although the ranges are subject to change 
given expected improvements in both data and model- 
ing, the ranges of parameters and correlations shown in 
Figure. [7] are remarkably close to expectations from less 
rigorous searches. 

The statistical procedures applied here represent a sig- 
nificant improvement to the state-of-the-art for compar- 
isons of data and models in the field of relativistic heavy 
ion physics. Previously, parameters were varied either 
individually, or in small groups. Figure [4] demonstrates 
the success of using emulators for this problem. Most 
importantly, the emulator techniques should scale well 
with increased data and increased number of parame- 



ters. Ultimately, the number of simultaneously varied 
parameters might increase to around 20, with the expec- 
tation that many, or perhaps most, of these parameters 
will not be significantly constrained by the data. Addi- 
tional parameters to which the model is insensitive does 
not increase the need for additional runs, as long as the 
parameters are varied in such a way that those parame- 
ters that are important are well sampled. It is expected 
that the additional parameters would decrease the effi- 
ciency with which the critical parameters are sampled, 
but that the additional number of model runs would not 
make the problem intractable. Adding more data should, 
hopefully, increase the number of principal components 
extracted from the model runs by providing additional 
discriminating power. Once the model calculations have 
been performed, the numerical cost of the statistical anal- 
ysis performed here was negligible, and adding more prin- 
cipal components should not cause any problems. Thus, 
the results of this study are promising, and encourage 
extending the scope to much larger data sets and more 
realistic models. 

Although the models used here represented the current 
state of the art, several improvements are necessary be- 
fore firm quantitative conclusions can be extracted. The 
following improvements require significant development, 
but are all tractable. 

• A flexible equation of state. For this study, the 
equation of state was fixed. There is some uncer- 
tainty involved with lattice calculations that should 
be accounted for with a variable equation of state. 
Additionally, it is of interest to address whether 
the equation of state is constrained by experiment 
alone, i.e. without relying on lattice calculations. 

• The bulk viscosity was set to zero here. Near T c , 
the system is undergoes a rapid change in micro- 
scopic structure, and the system may lose equilib- 
rium. As the system returns to equilibrium, en- 
tropy is generated. If the departure from equilib- 
rium is small, the effect can be accounted for by 
adding a bulk viscosity [17] [18] . If the departure 
is large, other approaches are possible, such as dy- 
namically solving for the mean fields [63] . 

• The chemical composition of the hadronic phase 
was set by the assumption of chemical equilibrium 
when the system reached a temperature of 170 
MeV. The chemical evolution can be improved by 
incorporating more inelastic process like baryon- 
antibaryon annihilation [64] [65]. Further, the as- 
sumption of perfect equilibrium at a fixed tem- 
perature should be relaxed by parameterizing non- 
equilibrium effects. 

• Although collisions produce many thousands of 
particles, the initial collisions involves only on the 
order of 100 nucleons. The finite number of original 
scattering centers leads to lumpy initial conditions, 
unlike the smooth initial conditions used here. If 
the model used here were improved to incorporate 
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FIG. 7. (color online) The distribution of acceptable values for each of the six model parameters are shown along the diagonal. 
The off-diagonal plots display the correlation between all pairs of observables. Four of the six parameters refer to the initial 
state (hence the "I.C." in their name) and the last two describe the shear viscosity. 
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FIG. 8. (color online) Twenty random points in parameter 
space were taken from the prior (upper panel) and posterior 
(lower panel) distributions. The temperature dependence of 
the viscosity to entropy ratios is clearly constrained by the 
statistical comparison with data, though the posterior distri- 
bution still covers a large variation. 



initial fluctuations, V2 would be more realistically 
modeled, and it would make it possible to consider 
fluctuations of the flow encoded in higher harmon- 
ics, i.e., ^3, V4 . . . [23] . 

Although the model used here can incorporate 
three-dimensional flow, for this study the calcula- 
tions were performed with the Bjorken ansatz. This 
approximation is reasonable for collisions at 200A 
GeV or higher [35], but full three-dimensional cal- 
culations are needed for lower energy, or for observ- 
ables away from mid-rapidity. Using low-energy or 
non-mid-rapidity data will also necessitate a more 
complex parameterization of the initial state. 

During the hadronic phase, pionic phase space be- 
comes highly filled at low p t . This affects spectra 
at the 10% level, which is neither crucial or non- 
negligible. Hadronic cascades can incorporate such 
effects by adding (1 + /) phase space enhancements 
to scatterings. This increases the numerical cost of 
the modeling at the factor of two level. 



Most of the improvements listed above would be accom- 
panied by an increase in the number of parameters. Vary- 
ing the equation of state will involve the addition of a few 
parameters. Since the bulk viscosity is not well deter- 
mined by lattice calculations, both it and its temperature 
dependence require parameterization. Non-equilibrium 
chemistry can be parameterized by adding fugacities for 
the initial state. For instance, the initial conditions away 
from mid-rapidity, or for collisions at lower energy, could 
necessitate a half dozen new parameters. At lower en- 
ergy, the dependence of the equation of state on baryon 
number is unknown and requires parameterization. The 
initial conditions at the LHC requires additional param- 
eters to encapsulate the beam energy dependence of the 
initial density and flow profiles. The lumpiness of the 
initial state involves setting a transverse size of the fluc- 
tuations. It is easy to imagine future analyses involving 
on the order of 20 parameters. 

The amount of available data for analyses such as these 
has swelled in the past few years. The beam-energy scan 
at RHIC provides all the observables analyzed here at a 
half dozen more energies. Additionally, Cu+Cu, Cu+Au 
and U+U collisions have been measured at RHIC. Addi- 
tionally, results from Pb+Pb collisions at the LHC have 
now been analyzed and published. Finally, higher flow 
harmonics, V2, . . . , v n , have also become available. 

Expanding the scope of the analysis to a larger range 
of beam energies and to include initial state fluctua- 
tions could increase the numerical cost of the calculations 
by two orders of magnitude. In the present study one 
processing core could perform a full model run for one 
point in parameter space in approximately one day. This 
would increase to being on the order of several weeks, 
or one month, if the beam energy scan, LHC data, and 
initial-state fluctuations were included. If the number 
of points sampled in parameter space were increased to 
a few thousand to better account for the larger number 
of parameters, the project would remain tractable, but 
would clearly require significant allocation of resources. 
The success and scalability of the methods presented here 
suggest that such an effort could help transform heavy 
ion physics into a more rigorously quantitative science. 
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